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ABSTRACT 


While solving internal flow equations near the 
inlet/starting stage, it is not always desirable to incorporate 
simplifying assumptions. In the past two decades some numerical 
schemes have been reported. in the literature, to solve 

multidimensional compressible Navier Stokes equations. A 

generalized implicit finite difference method with increased 
computational efficiency has been developed by W.R. Briley and H. 
McDonald [13. 

In this method the governing equations are replaced by its 
implicit time difference approximations. Terms involving non- 
linearities at the implicit time level are linearized by Taylor 
expansion about their solution, at the known time step. This 
results in a system of coupled equations which are linear for the 
dependent variables. Alternating Direction Implicit (ADI) 
technique is applied to the linearized equations. Here ADI 

technique developed by Douglas-Gunn has been applied. Application 
of ADI technique leads to equations which are one-dimensional. 
These equations are replaced by their finite difference 
approximations using the standard central difference formula. The 
finite difference equations when applied to successive grid points 
give rise to block-tridiagonal and tridiagonal systems of 
equations which are solved by standard methods. 

The method, illustrated above, is applied to solve the flow 
field through a rectangular duct being fed through a stagnant 
reservoir of fluid at one end and a specified pressure at the / 



downstream end. The primary aim is to predict the detailed 
compressible, three-dimensional flow field to establish the 
validity of the computer code, developed for this purpose. The 
computer code may be used for solving flow field in intakes and in 
solid rocket thrust chambers. Computations have been done by 
varying the input parameters to determine their effect on the flow 
field as well as the limitations of the method, if any. 
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1, INTRODUCTION 

1.1 General background 

A set of mathematical equations governing the flow of a 
compressible, viscous fluid was developed early in the 19^^century 
and subsequently became known as the Nav i er— Stokes equat i ons . These 
equations, in the most general form, include the' equations of 
conservation of mass .momentum and energy relating the pertinent 
dependent variables v i z . dens i ty , ve loc i ty . pressure and temperature 
which describe the flow of a compress i bl e . v i scous and heat 
conducting fluid, Since that time the basic form of the equations 
for describing the behavior of a continuous fluid has remained 
unchanged and form a set of non-linear coupled partial 
differential equat i ons . The work of subsequent research has been 
directed at finding approximate or analytic solutions for simple 
or special flow situations as well as to consider approximate or 
modified forms of these equations to predict turbulence .heat and 
mass transfer , shock waves etc.C8] 

There are basically three approaches or methods which can be used 
to predict the behavior of fluid motion. These methods may be 
grouped as <1) Experimental (2) Theoretical and <5) Numerical .The 
experimental approach has the capacity of producing the most 
realistic answers for many flow problems however the costs of 
measuring equipments etc. are becoming greater everyday as well as 
there is also the need to refine the measurement techniques for 
reliable measurements in complex flow cases Moreover in performing 
experiments it is not always possible to simulate the true 
operating conditions of the prototype. In the theoretical 
approach , s impl i fyi ng assumptions are usually made in order to rtraike 
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the problem tractabl e .The big advantage of the theoretical 

approach is that general information can be obtained, in many 

cases, from a simple express i on . Thi s approach is quite useful in 

preliminary design work since reasonable answers can be obtained 

in a minimum amount of time. However in many problems of practical 

interest analytical solutions are not possible. In the numerical 

approach, a limited number of assumptions are made and a high-speed 

digital computer is used to solve the resulting governing fluid 

dynamic equations. A parametric study can be easily carried out to 

study the variations in a given flow s i tuat i ons . C2] 

There are basically two numerical approaches , v i z . the finite 

difference method and the finite element method. [93In the former 

method, the governing differential equat i ons are approx i mated in 

terms of the dependent variables at a finite number of points in 

the flow field. There are many methods by which this can be 

achieved. In general, the partial derivatives of dependent variables 

at a point are replaced by expressions given in terms of their 

values at the neighboring locations and the distance between 

are solved 

them. This procedure gives rise to algebraic equations which. 

A 

simultaneously to obtain the values of dependent variable(s> at 
different locations. 

The finite element formulation is based on the integral statements 
of the conservation principles .The region under study is divided 
into a number of finite elements. The variation within the elements 
is generally taken as linear and the integral statements of the 
governing conservation equations postulate yields the integral 
equations that apply for each element .Minimization of the 
integrals is carried out in order to satisfy the conservation 



pr i nc i pi es , thereby obtaining the distribution of the variables in 
the reg i on . However in fluid flow problems finite difference method 
has been extensively developed and used during the past few 
decades . 

Problems of fluid mechanics may be divided into two 

categor i es . Steady state and unsteady state. The former relates to 
the circumstance when the variables at all points in the region of 
solution are independent of t ime . Unsteady state problems relate to 
the situation when variables in the region of interest vary with 
time, in addition to locat ion. Time-dependent solutions are of 
particular importance to manufacturing processes in which the 
thermal cycle undergone by the material is of i nterest . A1 so the 
start up and shut down of systems such as power 

plants , furnaces , ovens , rocket motors etc. invariably involve a 
consideration of transient conduct i on . Such systems frequently also 
have fluctuations in the boundary condi t i ons . Steady-state problems 
are of interest in the long range steady operation of systems, such 
as power systems , blast furnaces , ovens , reactors and storage 
un i ts . [ 91 

Finite difference schemes are further classified as Explicit 
scheme and Implicit scheme. An Explicit scheme is one in which only 
one unknown appears in the finite difference equation in a manner 
which permits its evaluation in terms of other known quant i t i es . An 
Implicit scheme is one in which more than one unknown appears in 
the difference equat ion .Th i s indicates that the formulation would 
require the simultaneous solution of several algebraic equations 
involving the unknowns. In both the Explicit and Implicit schemes, 
depending upon the step-s i ze , the solution may not always converge 
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Implicit schemes are more stable in the sense that they permit 
larger values of time step. Hence less computer time is required to 
obtain a so 1 ut i on .However Implicit formulation is more complicated 
to cast into a computer code.lt requires more calculations at each 
time level and hence requires larger computer time at each time 
step than that for explicit schemes. With larger values of time 
step, the truncation error will be more, and the use of implicit 
method to follow the exact transients may not be accurate . However 
for problems in which steady solution is required, the relative 
time-wise inaccuracy is not important. 

Several computer models for predicting three - dimensional, 
compressible flows have been reported in the literature during the 
last decade or so . Researchers at Los Alamos, USA, Imperial College, 
London, and other research workers have developed standard 
packages of algorithms such as TEACH code .MINT code ,etc., some 
of which are made available on commercial basis. 

1.2 Present, investigation 

In general, most of the flow field analyses are restricted to the 
steady state situations as well as for simpler geometries 
especially away from the inlet sect i on .However analysis of the 
internal flow field near inlet region is important for accurate 
prediction of the overall behavior of the fluid in several 
propulsion systems .During the burning of a solid propellant grain 
in a rocket thrust chamber, the combustion gases move downward at 
an increasingly higher speed before passing out of the nozzle. At 
the upstream end of the nozzle, the flow is developing essentially 
from stagnant conditions and the flow is fully choked at the 
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throat. The flow in the developing zone has low Reynolds number and 
also cannot be treated as a fully developed flow. For correct 
predi ct i ons , i t is necessary to deal with the flow field in this 
region which is of low speed and usually three -dimensional 
depending upon the grain geometry. For this purpose. a computer 
code for cold flow at low Reynolds number is developed to obtain 
steady-state solution of unsteady, three-dimensional, compressible 
Nav i er-Stoke ' s equation(s> based on a generalized implicit finite- 
difference scheme due to W.R. Briley and H. McDonald [1]. 



6 


2. PROBLEM FORMULATION 

2.1 Physical description of flow field 

IS 

The flow under cons i derat i on^ lami nar .subsonic .compessible in 
the entrance region of a straight duct having rectangular cross 
section. The duct is fed from a stagnant reservoir. The flow 
geometry has been shown in fig 2.1 

The flow consists of boundry layers on the walls and an inviscid 
core flow. At some downstream distance the boundry layers grow to 
fill the duct and the flow becomes fully viscous. 

The flow accelerates along the downstream and at some downstream 
distance the mean Mach no. becomes unity and the flow becomes 
choked .Depend! ng on the inlet Mach and Reynolds numbers .the flow 
may become choked either before or after, it becomes fully viscous. 
Here the duct length is chosen such that the flow at the exit is 
neither choked nor fully viscous. 

The flow is symmetric about the horizontal and vertical planes 
passing through the duct central ine .Therefore solutions have been 
computed only for one quadrant of the duct. Symmetry conditions 
have been imposed on the planes of symmetry. 

2.2 Governing equations s- 

The conservation equations for a compressible, viscous, heat- 

conducting, perfect gas, without body forces and external heat 

addition may be written in the following form. [33 
( Here i,j and l are being used as tensor indices ) 


Continuity equation 

a . a i 


at 

Momentum 


ax . 


& xi .y 

j 


equat i on 


a ( 0 u. y 
at ' 


ax 


= 0 


{ P U.Xl . 

«• J 


( 2 . 2 . 1 ) 


ap 

ax. 


\ 


+ 



T. . 

^ J 


< 2 . 2 . 2 ) 
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Energy equation 
a i 9 H ) ^ 


at 

where 


a i e M H '> _ dp 

ax . •’ ~ at 

j 


d < T. . 

^ ,, J 




(2.2.3) 


'• J 




2 J. a ■u^ 

— u 6. .— 3 — I 

9 V J 

err 

ax . 


r a XL, 

^ ( -a^r ^ 




Bx. 


J 

1 

. XI 

2 V V 


w = ^ + 

The pressure can be eliminated as a dependent variable 


means of the equation of state for a perfect gas 

p = 9 R r 


2.3 Botmdary conditions 
plane BC' G' F 

= 0 = V = ii> 

an 


dx 


= 0 


( no slip condition ) 

( adiabatic condition ) 


<2 


(2 

(2 


by 

2.4) 


3.1) 

3.2) 


dens i ty 
reduces 


can 

to 


be determined 

a ^ . a 


at 


from the 
9u 


dx 


= 0 


continuity equat i on , wh i ch 

(2.3.3) 


Plane A' D' H' E' 
a 9 _ Q _ du _ dv 
dx dx dx 


dw ^ dr 
dx dx 

( symmetry conditions) 


(2.3.4) 


Plane BA' E' F 

■u = 0 = V = ID (2.3.5) 


dT 


0 


(2.3.6) 


density can be determined from eqn. (2.2.1) 
a 9 ^ a 9^0 _ „ 

dt ay ~ ^ 


which reduces to 
(2.3.7) 


Plane c' D' H' G 

ay ay ~ ay ~ ay ay 


Plane A' BC' D' 

Since the duct is fed from a large stagnant reservoir it is 
reasonable to assume the flow upstream of the inlet plane 
as inviscid and adiabat ic .Thus constant stagnation 
temperature and pressure are assumed at the inlet plane. 
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T = T 

stag 


^■tgg 

P 


+ 


1 


2 C 

P 


2 


T 

• i €10 

_ 




-1 ) 


(2.3.9) 

(2.3.10) 


The velocity field at the inlet plane depends on 
sorrounding the inlet. — ^ '2 = 0 

has been used as a boundary condition which is 
with the velocity behavior observed experimental! 
boundary layers and inviscid core region. 

The other boundary conditions are 

du _ Q _ 

9s! ~ 3s 


conditions 
(2.3.11) 
compat i bl e 
y both in 

(2.3.12) 


Plane E' FG' H' 

The duct has an unobstructed exit into a large constant 
pressure reservoir. So it is reasonable to assume that static 
pressure at outlet remains constant 

i 


= g R T (2.7.15) 

and the other boundary conditions may be reasonably taken as 

(2.7.14) 




0 = 




^2 

a _ 
-_2 




2*4 Non di mens! onali sail on of equations 

The dimensional variables are normalized, with the reference 
quantities, in the following manner 

X = X/Lr ,y = y/Lr ,2 = ,U = ^t/Ur ,V = V^Ur ,W = W^Ur 

p = Q/p ,T = T/Tr ,t = tUr/Lr ,H = H/uf ,p = p/p U? 

r r 

This normalization leads to the following non dimensional 

parameters ; Mach no.. M = Ur/c ; Prandtl no., pr = ^^cp/k ; 
Reynolds no.. Re = p UrLr/p where c .the reference speed of 

r* 

2 

sound, is defined as c = j'RTr 

The fluid properties have been assumed to be constant. 

(Toverning eguatioiis 

Continuity equation 

3p 9 pu . 

9\ 9y. 

j 


0 


(2.4.1 ) 
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Momentum equation 

d pu. ^ d pu. u . 

V + V J 
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Energy equation 








a ( v.vT) 

dx . 


< 2.4 


g< ipT)_^ p u T>_ ^2, . 

at ax . Re pr 

j 


y<?'-l )M* 



p(V.^ 


( 2 . 
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( 2 . 4 . 18 ) 
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3. NUMERICAL METHOD 


3.1 Gonoral ourtlin© 

The governing equations are replaced by its implicit time 
difference approx imat i ons . Here backward difference scheme has been 
used for this approx imat i on . Terms involving non-linearities at the 
implicit time level are linearized by Taylor expansion about the 
solution at the known time step . This results in a system of 
coupled differential equations which are linear for the dependent 
var iabl es . 

Alternating direction impl i c i t < ADI ) technique is applied to the 
linearized equat i ons . There are many ways of applying ADI 
t echn i que . The ADI technique developed by Douglas-GunnC5] has been 
utilized here. The advantage with this ADI technique is that 
physical boundary conditions for the dependent variables can be 
use for the intermediate steps. ( Appendix ID. >[6,73 
The equations obtained are replaced by their finite difference 
approximations using the standard central difference formula. When 
applied to successive grid points and incorporating the finite 
difference equations for the physical boundary conditions 

block-tridiagonal and tridiagonal systems of simultaneous 
equations are obtained which can be solved by standard block 
elimination methods . [2 , 4 , 53 

3.2 Finit© diff©r©nc© ©quations for tb© govorning ©quations 

Continuity equation 

Implicit time difference approx imat i on of eqn. <2.4.1) at the 
unknown <n+l>*'^time step can be expressed as 


- 1 

dp 

dt 


+ 


d(pv ) 

dy 


1 


1 


^ . n4. 1 

^<pu> 

ax 


aipw") 

diz 


0 


<3.2.1> 
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Using backward time difference approximation .the above equation 
i s wr i 1 1 en as 


^ d(pu)"*i d(pv>'‘t‘ d(pw)’^*^ 

—Kt ^ W~ 5i— 


= ■ 0 


(5.2.2) 


It is convenient to work in terms of the difference of the 

variables at n^^ and (n+l>^^ time steps . So eqn. (3.2.2) is 
written as 


+ At + At 2^^ + At 


<3.2.3) 


p dx ay '■ a-z 

The linearization of eqn. (3.2.3) is done by Taylor expansion of 


<pu ) 
<pv) 
<pw) 


1 inear 

terms 

abo 

ut the solution at 

the 

known time 

step. That is 

n+ 1 _ 

(pu)"^ 

+ 

n n4- i 

p y>^ + 

n n«*> 1 
U V 

P 

+ 

o<At*) 



n+ 1 _ 

(pv)"^ 

+ 

n n+ 1 

P V'^ + 

n n-4* 1 

V 

+ 

o<At®) 


(3.2.4) 

n->* 1 _ 

<pw)^ 

+ 

n r>-»- i. 

P V'^ + 

r» n-»- 1 
W W 

p 

+ 

o<At^) 




Making use of approximations from eqns. <3.2.4) , eqn.<3.2.3> 

becomes 

n*l. , I L 9 , r, n+i , r> n+4 . , » j. ^ n+1 n n+1. . , 

Vf + At s—ip V + u v/ ) + At -=— <p + y yj ) + 

^p dx u p ay ^y ^p 


i . a . n n+l , n n+l 

At -s— <p y + WW ) = -At 

az w p 


f 3<pu) 
\ ax 


” + aipy ")^ ^ d<pw)'" 


dy 


az 


} 


<3-2.5) 

This gives a differential equation which is linear for the 
dependent variables .Now application of Douglas-Gunn ADI technique 
gives the following equations 


+ At 3 — <p V' + u V ) 
p ax ^ u p 


- At 


{ 


^(pu ) 


d(pv>” ^ aipyt'i'^ 

T 


ay 


az 

(3.2.6) 


} 
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** . , ^ 9 , r, *• , n •« . * 

W + At ■=— (p V + V W > = If/ 

V dy ^p 


<3.2.7> 


T'-'-l' . A ^ 9 , T> n+1 , n n+1. *• 

w + At -3— <p W + W V ) = W 

^p dz W ^p ^p 


<3.Z.8) 


Finite difference approximations for the above are 


- -Ta^ ^ * (- ]v. 


At u 


At p 


( 1 • .* r 

It-a^K^, = ^4-5^ - 




ay 


az 


(3.2.9) 


t j. 

At V . 


At v’' 


( ''j_, -j *« ^ . r "j*:! 1 •• A. f '"i-‘ 1 •• 

r 2 I 2 2iy J^p.^^ [- 2 Ay j^v._^ 


I 2 J "j.! 


V', 


(3.2.10) 


( At W, -4 

ill Iv,"** 

^ J p.-. 


+ V' 


n+ 1 


t At W, ^ j- lAi. f-7, -I 

iill + f iLli 

2 Az J ^k*i I ^ J '^k-i 


At p: 


( Ail -x 

■'‘« 1,,"“ 

2 i "k.. 


V', 


(3.2.11) 


X momentum equation 


Finite difference 
illustrated above 


equations 
with the 


are obtained by applying the 
exception that the term 


steps 
9i'7.V*> . 

sr-. IS 


evaluated at the known time step. This is done to 


avoid the ADI 
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treatment of mixed derivatives. 


The finite difference equations are 


At 

2 n 

- — T — ( U ) , 

2 Ax v -1 


At T 


2 Ax y 


T1 

—*2 ]v/* + u” V'* + f 


At 

■ 5 — T— <U > . + 

2 Ax + 1 


At t: 


2 Ax 


H... * (■ 


2 Til r> . , 

p. u. At 

V - i V - ± 

p— 


Re Ax 




2 At 
Re Ax" 


2 u'' At At 


Jv'u, I Z Ax -- 


■j * , r ) * 

r 2 Ax ^ 


Re Ax 


At p. 


2 Ax y 


i + l I « 

J^t. 

^ ^ v + dL 


aipu^^^ 
ax ^ 


aipus/y. ^(puw>. 


_L.2 ^ 

;k ax 


1 r ^2 n 


n 1 dCV 

i 5 ax 


) } 


(5.2.12) 


n n 

At U . V . 


. , n n 

At U . V . 


2 Ay 


v-v ^AAL UV’V 

yp^., ""pi [ 2 Ay 


. T> n 

At p . V . 

J-i J-i 

2 Ay 


Re Ay 


1 «» 

J^u. 

J j-i 


n 2 At 

> . + —5 r 2 

j Re Ay 


1 »* 


. . n n 

At p Y 

j-i-i 

2 Ay 


At p^ 


ce Ay 


> *• X f 1 ** 

J^u . 2 Ay j^v . 


. . n n 

At p . u . 


( AJt L . U . ^ ^ 

Iv/** 


n * , n • 

P . F + ^ 

J y. j p. 

j j 


<5.2.13) 
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At u, w, 
k - r 


( 4U t U, W, -v 

k~i k "l I n4.i 

- ___ 


+ u" + 


At u, w, 
k - 


... 


(■ 


r> n 

At p, W, 

•^k-i k - 1 

Z~Ez 


At 


Re Az 


K:.^ [ 


n ^ 2 At 

^k Re Az* 


1 n^± 

K 


f 

n n 

Pk + i'^k+i 

At -x 

2 + 

f_ 

At 

n n 

1 2 

Az 

Re Az J^u 

J k + l 

1 

2 

Az J 


f At p, u, -» 

1 k + i k-*-l I n+4 

I 2 Az 

k J k + i 


n *» , n »» 

V'.. + k.,_ V', 


k " U. 




(3.2.14) 


y momentum equation 

Finite difference equations are obtained by applying the steps, 

a(V.v^) . 


illustrated above. Here the term 


dy 


is evaluated at the known 


time step. 

The finite difference equations are 


n n 
At u. V. 

V * " 


( XiL U. V. "V ^ 

2 J Pi-. 


+ V . + 


. . n n 
At U. V 


I — r -^ — Jv * 


(- 


. . n n 
At p. V. 

2 Ax 


. - n n 

P. . . V, 


* . n n 

At p. u. 


]yr + f ly* + f- 

I^U I 2 Ax Pu. I 2 Ax 

-'V-i -' 1 . 4-1 V 


At 


Re Ax' 




n ^ 2 At 

fte Ax* 


. . n TV 

At p. u. 


•v _ r ■“T. u. 

K ^ 


1 • _ A* f d<puv)" a<pv*)*' aif 

Re Ax* j^'^v + 4 V 
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1 + 
rV sr ‘ * 




< 3 . 2 . 15 ) 


2 Ay 


/ 2 . r* 

<V > . 


At T” 

j_ 

Z Ay r 


n 

j_i l*» . n »* r 

^>.2 |V + V . V +1 

W J ^j -1 J I 


At 

, 2 n 

■= — r— (V > , + 

2 Ay j+i 



f 

2 p V . At 

j-i 1-i 


f " . 

Vp. 

1- 

2 Ay 

Re Ay= l''v^ * 

j 4-1 



2 At 


2 p'' u" At At 


]<: * ( 


A J. ^ 

At p . 


Re Ay 


1 /* . f- ... r . ri - s 1 /* 

J^v . I 2 Ay y M Pt . 

J j+i ^ J j-i 


r* 

At p; 


+1 I Ay Y 


j*i 1 *» 

rr. 

J j*± 


n • , n • 

P . V + V . V 
""j %. J V. 


( 3 . 2 . 16 ) 


^ n*l n n + 

— \w + V, y 

FAi J ^k-i ^ 


. , n n 

>. ^ At V, W, X 

I r>-»-l , n , | k *4- 1 k-*.il r>*^ 1 

\ * i ' K.. 


F~a 5 


Re Az 


1 T>+ 1 

jV 

7 k - 1 


^ r ^ u. 2 At 

+ ^k + -RTAz^ 


1 n-i- i 

J’'Vk 


^ ^ * X 

M p w. At 

^k •*■ i. k J. 


T’^A^ 


Re Az 


1 n 4 .i 

J k-*. 1 


. . n n 

At p, V, 
•^k-l k 


2 Az 


1 k-l 1 n+i 
— 

J k-l 


At Pu*±''k*± 

Faz 


1 n + 1 

K. 


n , r. ** 

p, y + V, y 
k k 


( 3 . 2 . 17 ) 


2 momentum equation 
d(V ) 

Here the term ^ — is evaluated at the known time step. The 


I A> A 


Hi -F-f :prpnr!e eauat i ons are 
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. . n n 
At U. W. 


( iUk L u . w. ^ ^ 

t.Zl- ]y* 

2 J 


, n * 

+ W. + 

P. 

1 . 


AX 

At u. w. 


L U . W. A _ 

ili—L** lu/ 

^ J ^i.it 


. , n n 

At p. W 


AX ^ ^ 

^ XkL p. w. . "I ^ ^ At p. w. V ^ HZ p, u, 

I v-i t-i I * I '^V4.iv + 1 }♦ ^ I *^v-li-l 

I 2 Ax I 2 Ax r 2 Ax 


* . n 1 ^ 

At p. u. 


At 


Re Ax 




n Z At 

+ 2 


A . n n 

At p. U. 


Re Ax 


Jv I — — 


At 


Re Ax 


jv,* = At [- 

-I X X £ t 


. . , d<puw>, d(pvw>. 3{pw >. 

At •{ - s — V - -s V - -s — V 

dx dy dz 


1 a<pT) 1 

y M* dz '■ Re 


(r»2 ^ 

V w. 


1 d(^.V 

? dz 


>rj} 


(3.2.18) 


A X n 

At V . W . 


( XX U V . VW . ^ ^ 


+ + 

a V. 


AX ^ ^ 

At V . w 


I 2 Ay 


A X ^ 

At p . w . 


. . n n 

At p . w . 


. , ri n 
At p . V 


( XIT. p . W . ^ -V ^ XX L p . W . -V f XX L p . V . 

- J-* Iv,** + I Iv/ + I 

2 J '^j-i I ^ J ''j*i I ^ 


At 


Re Ay 


1 ♦ilt 


r " 


Z At 1 

^ ( 

J^w 

J J-1 

+ 


+ 

Re Ay* J^w. 

" 1 


. n r> 

( At p . V . 
2 Ay 


At 


Re Ay ! ' w . 

* jfi 




n • . n * 

P . W + W . y/ 

^ '^j ■’ 


(5.2.19) 


(- 


At 

T-a^ 


At T 


k-l 


2 Az r M 


1 r>*i 

K-, 


, n n* 1 
+ W, W 

‘ '’i 


( 


At 

T-T- + 

2 Az k + i 
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At T 


2 Az 


*1 "J n*! f 

K.. * [- 


2 ?*^ Tip . , 

P, w. At 

•^k-l k-i 

T~Az 


At 


Re Az 


)s:. ^ ( 


* 


2 At 
Re Az' 


2pr. w^.'At At 


r •- H, 

J r.-1-l. 1 ^k-t-l k+1 

JV ^ I — A i- - 


•V ^ A t f> -v 

I n*»-l , I k-i I r» 4- 1 

r 2 Az r m" 


Re Az 


( All x 

k 4-1 I n-»- 1 

2-ArFM" 


n *• . n •• 

= p y/ + w, V' 

•^k ^w, k ^p, 

k k 


(3.2.20) 


Energy equation 

While obtaining the finite difference equations for the energy 
equation, the diffusion term *p is evaluated at the known time step. 

This is done to avoid , the finite difference equations , from 
getting much complicated. 

f- u" T” 1/ + f t" + AK^'-Df Ji-VT" 1/ + 

2 Ax v-i v-ij^p^ ^ L*' '-dx-'vv j 


f At n _n 1 * , f 

[ TAJ * [- 


At n _Tr> 

P: . T, 


2 Ax "^i-l i-i 


At(r-1> n _n 1 * 

2 Ax Pi ^i J'"u. 


( At n — n 

2 Ax ^v*i v*i 


AtCyl ) 
2 Ax 


pr^rR ^ (- 

V 4. 1 ^ 


At 


n n 

P: . 


2 Ax i-l 


At ^ 
fte pr Ax'* 




. z r 1 * 

fte pr Ax* 


( 


At n n 

P: . . 'J: 


y At 


2 Ax + i i + i Re pr Ax 


jr* = At 

%.4*1 V 


a(puT ) ^ a(pvT > 

^x ** ^y ^ 
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^(pwT ) . 
_ 


(y'-Dp" <V.V^)" + 

V X, \ 


^ (V^T)"' + 

Re pr V 


.n 

Ri 


} 


0.2.21) 


r At n _n 1 »* J 


3\/ 

? _ri 1 •» 

1" 2 Ay ''i-.^-.J'"p-.. ^ 1 

1 T^ + At<j'-1)( 

ay 

-1 T Iv/ + 

Jj J JV. 

f At n _n 1 •* . J 

1 V T 1*1/;/ + 1 

At n 


At(?'-1> n 1 

[ 2 Ay j + ^ 1 

2 Ay Pj-t' S-± 


2 Ay J 


( 


j. t t'^ j. AKj'-l) 

* ' Ta9 ^ - 2 aV ^ 


^ I- 


At n n 
2 Ay ^j-i''j-i 


y At 

Re pr Ay" 


]'"**-, *('’j + 


Re pr Ay 


1 mm 


At r> n 

P: . 


2^ At 


2 Ay j+i j+i Re pr Ay 


j^T. 

J + i 


^ m ^ <ak. M V 

= V'y + V'p 0.2.22) 

j 


( At n _n 7 n+i T n , ... i n r ^ 7 , 

■ IS? "k- I !> C 5i-3k\ J'"p^ + 


r At n _n 1 n+l 

x 

f- t" 

Att^^-l) 


[ 2 Az k*± k-RiJ^Pj^^^ 

1 

[ 2Az^k-i‘k-i 

2 Az 

^k 'k J 


^ r At n _n ^ At(>'-1) 

+ • P ,. . T,_ _ + ^ 


2 Az k*£ k + i 2 Az 


n _n 1 n+± . f 

'’k \ jv * (■ 

k-i-t k 


At 


n n 

Pu 


2 Az '^k-l k-i 


r At 

1 3 

Re pr Az 


In+l , f n ^ . 2 ^ At 1 n-4-1 

J’'\-. 1'°’' Re pr A z^ J'"t^ 
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( 


At n n 
O W 

2 Az k+i 


^ At 

Re pr Az^ 



i 

-»> 1 



mm 



+ 





(3.2.25) 


3*3 Find to difforonco equations for the boundary conditions 

For obtaining the finite difference equations for the boundary 
conditions, the steps, which have been already followed for the 
governing equations, are applied . Three point, one sided , 

difference approximations are used for the boundary conditions. 


Plane B C' G' F 

Finite difference approximations for eqns.(2.4.6) are 


n + 1 « n -*- 1 n -»- 1 

w = 0 = V = V' 

V w 

i=0 O O 




Finite difference approximation for the equation (2.4.7) is 

3 V' j - 4 yj = 0 (3.5.2) 

i =o i 2 

Density is determined from the equation (2.4.8). Implicit time 
difference approximation of this equation is 


dp ^ a(pu ) 


at 


dx 


. .. d(pu) 

W + At 
P 


n *- 1 


9x 


= 0 


Linearization followed by space differencing of the above equation 
gives 

3 At n ] , f 4 At , n ] ^4-1 f At n 

L ■ 2 Ax ^i=o [ 2 Ax i yp^ [ 2 Ax 2 

I’ r-E? '’o Jv t r-E? J'"u, ^ I TTz '’z 
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.. d(pu) 
At -s-t— o 
3x 


<5 


Plane A' D' H' E' 


Finite difference approximations for the equation (2.4.9> are 


n-i- 1 

^ V =L-2 

- 4 

r»*«. 1 

+ 

3 

n 4. 1 

n+ 1 

W 

u. 

V =JL-2 

- 4 

n-4- 1 

U 

L-l 

+ 

3 

n+ 1 

¥ 

U 

L. 

n-»- 1 

V 

i =L-2 

- 4 

n4. 1 

¥ 

V 

L- 1 

+ 

3 

n4- 1 

W 

V 

L 

n*»- 1 

¥ 

\ =l--2 

- 4 

n+l 

L- 1 

+ 

3 

r»+ 1 

W 

W 

L 

r>4- 1 

^T. 

- 4 

n4. 1 

+ 

. 3 

1 

I. 


p + Ap - -Jp 

JL - 2 L - 1 


u + 4u - 5u 

3L-2 L-1 I 


V + 4v - 3v 

1--2 L-l 1 


w + 4w - 3w 

L-1 L 


- T + 4T - 3T 

L-2 L-l L 


(5 


(3 


(3 


<3 


(3 


Plane B A' E' F 

Finite difference approximations for eqns . ( 2 . 4 . 10 > are 


n^-j. _ n-*-! n+1 

W = Q zz w = 

j = o o o 


? V 


n -»- 1 


jsrO 


n + 1 

^ V^j + 

1 


n-i- 1 

2 


(?, 


Finite difference approximation for the equation (2.4.11) is 


(3.1 


Density is determined from the equation (2 . 4 . 1 2 > . Imp 1 i c i t 
difference for the equation (2.4.12) is 

. n+ i 


r>+ 1 n 

p - p 

At 


Sips/ ) 

w 


= 0 


n+ 1 

y/ + 

P 


At 


_ , . n* i 

S i py ) 

dy 


0 


.5.3) 

.3.4) 

.3.5) 

,3.6) 

.3.7) 

.3.8) 

3.9) 

k 10) 

t ime 


Linearization followed by space differencing of the above equation 
gives 



2 ? 


fl - ^ 

At n 

1 n-i-l 

f ^ 

At n 1 n-^ 1 


r At n 1 

V ^ 

— T V . 

Ay j=o 

Vp 

[t 

< 

+ 

1 2 Ay ''2 j 


1 n+ 1 

K 


f 3 

At n 1 n-»-i 


r ^ 

At n 7 1^+1 


f At n 1 

1" 2 

> 

< 

0 


u 

a 7 '"i 


[ 2 Ay ^2 J 




n+ 1 
2 


. . a(pv) 
At -K— — o 


ay 

Plane C' G' H' D' 

Finite difference approximations for the equation (2.4.1?) are 


(?.?.11> 


n**- ± 

'^4=M-2 

- 4 

n + 1 

M-1 

+ 

« n + 1 
% 

= 

n 

- ^M-2 

+ 

X 

_ n 

- 

(3.3.12) 

n+ J. 

W 

U . 

j=:M-2 

- 4 

n+1 

M- 1 

+ 

- n+1 

3 

U 

M 

= 

n 

“ U 

M-2 

+ 

40*" 

M- 1 

- n 

- 

(3.3.13) 

n4 . 1 

W 

V . 

jsM-2 

- 4 

n + 1 

W 

V 

1 

+ 

- n + 1 

3 V' 

V 

M 

= 

n 

- V 

M-2 

+ 

X ^ 

4v 

M- 1 

- 3v"‘ 

M 

(3.5.14) 

n+1 

¥ 

w . 

j=M-2 

- 4 

n + 1 

¥ 

1 

+ 

1 n + 1 

3 V' 

^w 

M 

= 

n 

- w 

M-2 


4w"' 

M- 1 

- 3w’^ 

M 

(3.3.15) 

n-i“l 

^T. 

J=M-2 

- 4 

n + 1 

M- 1 

f 

- n + 1 

5 

M 

= 

- 

M-2 

+ 

4T" 

M-1 

- 

(3.3.16) 
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From equation (2.4.15) 

P 

e i oirg 


n+1 _n + 1 / 

^ 1 

- .j.n+1 -j7'/(j'-l) 

n _n 

p T 

f r 1 

2^ M 1 

T 1 

2^ M 

• J—— 1 

• V 

• 1 


• i a^f 




Linearization followed by finite difference approximation of the 
above equation gives 
1 


n-i* at 


r» 

P 

k = 0 n-H 

V't 


(5.5.17) 


k = 0 


ksO 


From equation (2.4.15) 
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T = T"*‘ + i r- P («=>"** = T" + ‘ ; .u 

stag £ Z 


/ 2 . n 

(W ) 


Linearization followed by finite difference approximation of the 
above equation gives 


. -...2 n n^-i 

^ k=0 ^W, 


+ y 


n*«- 1 


<?. 3. 18) 


k=o 


k = 0 


Finite difference approximation for the equation (2.4.16) is 


V'. 


r>-*- 1 


W 


- 2 


k = 0 


n+ 1 
W 


n-»- 1 
W 


V' 


W + 2w - W 

O i 2 


(3.3.19) 


Finite difference approximations for the equation (2.4.17) are 


3 V' 


n+ 1 


U 


- 4 


n-*“ 1 




n4. j. 


*T ” t A ^ ' 

= -3u + 4u ~ u 


(3.3.20) 


k = o 


3 


r>+ 1 


- 4 y; 


n-*- i 


4- = -3v" + 4v' 


(3.3.21) 


k = 0 
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From equation (2.4.18) 

. n+ 1 


• i at i c 


(p ry 

r M 


<p Tj 
r M 


Linearization followed by finite difference approximation of the 
above equation gives 


_n n-*- 1 

^k=N 


n n-«- 1 

P V'j 

M 


(3.3.22) 


Finite difference approximations for the equation are 


V', 


_ n+ £ 
U 


2 y/, 


k = N-2 


n-i" 1 


4 y/ 


n-i* £ 


N 


r> , ^ n n 

- u f -Zu - u 

H-2 H-£ N 


(3.3.25) 
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n**- i 
w 

v 

k = N-2 

2 

n-i- 1 

W 

V 

K- t 


n-i* 1 

V 

N 

( 

1 

II 

+ 

2v 

N- ± 

n 

- V 

N 

<5. 3. 24) 

n+ 1 
y/ 

ksH-2 

2 

n-*- 1 

V' 

w 

N~ 1 


1 

y/ 

W 

N 

n 

= “ W 

N-2 

+ 

2w^ 

N- 1 

n 

- VNi 

N 

(3.3.25) 

n+ 1 
vj 

k = N-2 

2 

n+l 

N- t 

+ 

n-*- 1 

N 
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3# 4 Solution proceduro for the finito difforenco equations 

Careful study of the difference equations reveals that, for each 
ADI step, it is required to solve one block-tr i diagonal system of 
equations and two simple tridiagonal systems of equations. 

During the first step of the ADI procedure (ie. marching along x 
axis), only the derivatives with respect to x and t appear in the 
implicit difference equations. In the continuity, x momentum and 
energy equations, these implicitly treated terms contain p ,u and 

T,but not V and w. Therefore the difference equations from these 

« » m 

three equations can be solved for p , u , and T .Having obtained 

« at at 

the values of p ,u , and T the difference equations for y and z 

<4(' 

momentum equations can be solved independently for v and w since 
y and z momentum equations are uncoupled with respect to v and 

■Hi 

w .Solution of V requires the solution of simple tridiagonal 
system of equat i ons . S imi lar 1 y solution of w also requires the 
solution of simple tridiagonal system of equat i ons .Th i s special 
nature of coupling improves the computational efficiency. 

During the second ADI step (ie. marching along y axis) the 
difference equations from the continuity , y momentum and energy 
equations are solved as coupled equations for p .v ,and T 

m# IK# 

Having obtained the values of p ,v ,and T the difference 
equations for x and z momentum equations can be solved 



26 


iK0 iimi 

i ndependen t 1 y for u and w since x and z momentum equations 
are uncoupled with respect to u and w 

During the third ADI step <ie. marching along z axis) the 
difference equations from continuity ,z momentum and energy 
equations are solved as coupled equations for p"** .and 

and difference equations for x and y momentum equations can be 
solved independently for u”**and v"**. 

First AD I step 

From equations (5-2.9) , (5.2.12) ,(5.2.21) and (5.5.1) ,(5.5.2) . 

(5-5-5) the following system of of equations is obtained 
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The system of equations (3.4.1) gives the values of p ,u and T at 

all the grid points except for j = 0 , M and k = 0 , N , for which 
the values are calculated directly from the boundary conditions. 
Equations < 3 . 2 . 1 5 ) , < 3 . 3 . 1 > and (3.3.7) give a system of equations 
represented as 
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System of equations <3«4,2> gives the values of v at all the grid 
points except for j = 0 , M and k = 0 , N for which the values are 
calculated directly from the boundry conditions. 

Similarly equations (3-2.18) , (3.3.1) and (3.3.7) give a system 
of equat i ons , s imi lar to ( 3 . 4 . 2 ) . Appl i cat i on of boundary conditions 
gives the values of w at all the grid points. 


Second ADI step 

Here equations (5.2.10), <5.2.16), (5.2.22) along with equations 
(5.5.9), (5.5.10), (5.5.11), (5.5.12), (5.5.14), (5.5.16) gives a 
system of equations of the form (5.4.1). 

Equation (5.2.15) along with equations <5.5.9) and (5.5.15) gives 
a system of equations of the form (5.4.2). 

Equation (5.2.19) along with equations (5.5.9) and (5.5.15) gives 

a system of equations of the form (5.4.2). 

«« mm •• .„•* . .... 

Thus the values of p ,u .v .w ,T are calculated at all the 
grid points except i = 0 , L and k = 0 , N , for which the values 
are calculated directly from the boundary conditions. 



?1 

Th i rd AD I step 

Here equations (5.2.11), (3.2.20), (3.2.23) along with equations 

(3.3.17), (3.3.18), (3.3.19), (3.3.22), (3.3.25), (3.3.26) gives a 

system of equations of the form (3.4.1). 

Equation (3.2.14) along with equations (3.3.20) and (3.3.23) gives 
a system of equations of the form (3.4.2). 

Equation (3.2.17) al ong with equations (3.3.21) and (3.5.24) gives 
a system of equations of the form (3.4.2). 

Thus the values of . u"^*^ . v’^* * .w^** .I’^'^^are calculated at all 

the grid points except i = 0 , L and j = 0 , M , for which the 
values are calculated directly from the boundary conditions. 

3« 5 Tridiagonal form 

Systems of equations (3.4.1) and (3.4.2) may be easily reduced to 
block-tr i d iagonal and tridiagonal systems of equations 
respect i ve ly .This is done so that standard algorithms may be 
applied to solve the system of equat i ons . Th i s increases the 
computational efficiency. ( Appendix I and Appendix II ) 


System of equation (3.4.1) is reduced to block-tr idiagonal form if 
the matrices E and F can be eliminated. E is eliminated in the 
following manner. 
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multiplying equat i on (3 . 4 . 3 b) by E C”* and then subtracting it from 
equation gives 
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are transformed to 

eliminate F . 


System of equations (5.4.2) is reduced to tridiagonal form if e 
and f are eliminated.e is eliminated in the following manner. 
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then subtracting it from (5.4.4b) 
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Therefore if e has to be eliminated then d , a and r are 

o o o 

transformed as follows 



where ^ and r"^ are the transformed values of d^ . a^ and r 

o o o o o o 

respectively. 


3. ^ Init-ial condit-ions 

Uniform flow i n the axial direction has been assumed as the 
initial condition. No slip conditons have been applied, where 



app i cab 1 e . L i near variation of density along axial direction has 
been taken. 

Stagnation pressure and stagnation temperature are specified and 
downstream static pressure is impulsively reduced to its specified 
value. The duct geometry has been taken as x^= 1 ” 


0.5 . 
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5. RESULT AND DISCUSSIONS 

le numerical scheme, as described in previous chapter, has been 
iplemented, with the help of computer, to solve the flow field 
irough a rectangular duct, fed through a stagnant reservoir of 
uid at one end and a specified pressure at the downstream end. 
le primary aim is to predict the detailed compressible, three- 
mensional flow field to establish the validity of the computer 
ide developed. The computer code developed may be used for 
riving flow fields in intakes and in solid rocket thrust chambers 
rmputations have been done by varying the input parameters to 
itermine their effect on the flow field as well as the 

mitations of the method, if any. 
feet, of different parameters 

gures 4.1 and 4.2 show the variation of axial velocity(w>, with 
at a fixed value of y at the upstream and downstream ends for a 
irtain values of input parameters viz. p , T , p , Re and M 

•log atag • 

lese curves indicate the typical characteristics of boundary 

tyer profile. Also there occurs thickening of the boundary layers 
» the downstream direction. This suggests that the flow is 

ive lopi ng. 

»e variations of centreline axial velocity and pressure along 
cial direction are shown in figures 4.3 and 4.4 respectively. The 
iow is found to be accelerating. 

1 figures 4.5 and 4.6 variations of centreline velocity and 

intreline pressure respectively have been shown for a higher 
ilue of p^. Also centreline velocity is less for higher p^. 
gures 4.7, 4.8 and figures 4.9, 4.10 show the effect of p^^^ on 


ie f 1 ow field. 


3S 


igures 4.11, 4.12, 4.13, 4.14 and figures 4.15, 4.16, 4.17, 4.18, 
ihow the effect of Re on the flow field. It is observed that 
ligher the Re, lesser is the boundary layer thickness, 
igures 4.19, 4.20 and figures 4.21, 4.22 show the effect of Mach 
lumber on the flow field. 

Effect, of mesh size 

lifferent mesh sizes were used to compute the results. It is seen 
;hat the solutions are in satisfactory agreement. However the 
:omputer time increases considerably for finer mesh sizes. system 
:ime increases from around 4 seconds for 6x6x6 grid to around 30 
seconds for 6x6x21 grid. The system times are for ’tsrx*. Here it 
s noteworthy that system times are higher for higher values of 
(e . 
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5. FUTURE WORK 

Deen mentioned that, computer time with finer mesh sizes. 
Ter large as compared to that with coarse mesh sizes. Also 
ocity gradients near the solid boundary are very large 
larly for high Re. These suggest that it is desirable to 
ce variable mesh size grid such that there are more number 
points near the solid boundar i es . Th i s may help in cutting 
e computer time. 

since in most real situations, the flow field is 
int , an appropriate turbulent model may be incorporated in 
jverning equations and the computer code modified 
i ngly . 
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Appendix I 

Solution of Tridiagonal system of equations 



r the solution of above tridiagonal system equations 
gorithm " is given as 


r i = 1 to L 


d. i - 1 

V - i 
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a* 9 . 

i 4.1 
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Thomas 



Appendix II 
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Solution procedure for block-tri di agonal system of equations 

C4.5] 



'e IB. A. C are <n y n) matrices, R is a vector <n y 1 ), ^nd G 

the column vector (n y 1) of variables. Solution procedure for 
! above system of equations is the following. 
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APPENDIX III 

Alternating direction In^licit Technique 


A = { F } 


(3.1) 


Here C is the column vector of dependent variables. £ , & , & 

X y * 

and A are the row vectors, ft , consist of partial 

x y z 

differential operators with respect to x, y and z respectively. 

ADI technique developed by Douglas-Gunn [1,6,7] when applied to 
the above partial differential equation gives 
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